load necessary packages

library(tmap)
library(mapview)
library(raster)
library(rgdal)
library(viridisLite)
library(RColorBrewer)

set working directory

path <- "/home/ede/tappelhans/uni/marburg/lehre/2016/ss/PS_global_change/code"
setwd(path)

read and prepare data

lyr <- ogrListLayers("../data/gadm_deu/gadm_deu.shp")
gadm <- readOGR("../data/gadm_deu/gadm_deu.shp", layer = lyr)
## OGR data source with driver: ESRI Shapefile 
## Source: "../data/gadm_deu/gadm_deu.shp", layer: "gadm_deu"
## with 403 features
## It has 15 fields
fls <- list.files("../data/tmean_deu", pattern = glob2rx("*.tif"),
                  full.names = TRUE)
stck <- stack(fls)

tmean <- mean(stck)/10

tmean_gadm <- extract(tmean, gadm, fun = mean, sp = TRUE)

quick look with mapview

mapview(tmean_gadm, zcol = "layer", alpha.regions = 1, legend = TRUE)

quick plotting with spplot

spplot(tmean_gadm, zcol = "layer")

using tmap package to produce cartography-grade maps in R

tm_shape(tmean_gadm, unit = "deg", unit.size = 1) +
  tm_polygons()

tm_shape(tmean_gadm, unit = "deg", unit.size = 1) +
  tm_polygons(col = "layer")

pal <- colorRampPalette(rev(brewer.pal(9, "Spectral")))

tm_shape(tmean_gadm, unit = "deg", unit.size = 1) +
  tm_polygons(col = "layer", palette = pal(6),
              auto.palette.mapping = FALSE)

tm_shape(tmean_gadm, unit = "deg", unit.size = 1) +
  tm_polygons(col = "layer", palette = inferno(9),
              auto.palette.mapping = FALSE,
              title = "Mean annual \ntemperature [°C]",
              border.col = "grey20", border.alpha = 0.3,
              breaks = seq(5, 11, 0.5)) +
  tm_legend(legend.frame = TRUE, legend.bg.color = "grey") +
  tm_style_grey(bg.color = "grey90", frame.lwd = 2) +
  tm_compass(position = c(.01, .055), color.light = "white", size = 3) +
  tm_scale_bar(position=c(0.03, 0.015)) +
  tm_layout(legend.title.size = 0.8)

change margins to avoid overlapping legend

tm_shape(tmean_gadm, unit = "deg", unit.size = 1) +
  tm_polygons(col = "layer", palette = inferno(9),
              auto.palette.mapping = FALSE,
              title = "Mean annual \ntemperature [°C]",
              border.col = "grey20", border.alpha = 0.3,
              breaks = seq(5, 11, 0.5)) +
  tm_legend(legend.frame = TRUE, legend.bg.color = "grey") +
  tm_style_grey(bg.color = "grey90", frame.lwd = 2) +
  tm_grid(projection = "longlat", labels.size = 0.5, lwd = .25,
          labels.inside.frame = TRUE) +
  tm_compass(position = c(.01, .055), color.light = "white", size = 3) +
  #tm_scale_bar(position=c(0.03, 0.015)) +
  tm_layout(legend.title.size = 0.8,
            inner.margins = c(0.02, 0.15, 0.02, 0.02)) +
  tm_xlab(text = "Longitude") +
  tm_ylab(text = "Latitude")

World example

data("World")

tm_shape(World, projection = "robin") +
  tm_polygons(col = "HPI",
              auto.palette.mapping = FALSE,
              palette = "RdYlGn",
              title = "Happy Planet Index") +
  tm_format_World(inner.margins = .02, frame = FALSE) +
  tm_style_classic(earth.boundary = c(-180, 180, -87, 87))  +
  tm_legend(position = c("left", "bottom"), bg.color = "grey95", frame = TRUE) +
  tm_credits(c("", "Robinson projection"), position = c("RIGHT", "BOTTOM")) +
  tm_grid(projection = "longlat", labels.size = .6, lwd = .25)

tm_shape(World, projection = "robin") +
  tm_polygons(c("HPI", "gdp_cap_est"), auto.palette.mapping = c(FALSE, TRUE),
              palette = list("RdYlGn", "Purples"),
              style = c("pretty", "fixed"), n = 7,
              breaks = list(NULL, c(0, 500, 2000, 5000, 10000, 25000, 50000, Inf)),
              title = c("Happy Planet Index", "GDP per capita")) +
  tm_format_World(inner.margins = .02, frame = FALSE) +
  tm_style_natural(earth.boundary = c(-180, 180, -87, 87))  +
  tm_legend(position = c("left", "bottom"), bg.color = "grey95", frame = TRUE) +
  tm_credits(c("", "Robinson projection"), position = c("RIGHT", "BOTTOM"))

See this online tutorial for more possibilities to plot spatial data in R.